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Owing to their favorable scaling with dimensionality, Monte Carlo (MC) methods have become the 
tool of choice for numerical integration across the quantitative sciences. Almost invariably, efficient 
MC integration schemes are strictly designed to compute ratios of integrals, their efficiency being 
intimately tied to the degree of overlap between the given integrands. Consequently, substantial 
user insight is required prior to the use of such methods, either to mitigate the oft-encountered 
lack of overlap in ratio computations, or to find closely related integrands of known quadrature in 
absolute integral estimation. Here a simple physical idea — measuring the volume of a container by 
filling it up with an ideal gas — is exploited to design a new class of MC integration schemes that 
can yield efficient, absolute integral estimates for a broad class of integrands with simple transition 
matrices as input. The methods are particularly useful in cases where existing (importance sampling) 
strategies are most demanding, namely when the integrands are concentrated in relatively small and 
unknown regions of configuration space (e.g. physical systems in ordered/low-temperature phases). 
Examples ranging from a volume with infinite support to the partition function of the 2D Ising 
model are provided to illustrate the application and scope of the methods. 



I. INTRODUCTION 



To introduce and place the ideas of the present contribution in context, consider the paradigmatic problem of 
estimating the volume V of a given region as shown in Figure [T] In its most rudimentary form, MC estimates of 
V proceed using the "hit-or-miss" idea, whereby the user designs a reference region of known volume Vq that fully 
overlaps with V, and draws random points uniformly distributed in Vq (cf. Fig. [I] top panel). An estimate of V 
can then be obtained from V/Vq — f, where / is the fraction of such points that fall inside V. As is immediately 
apparent, the efficiency of this simple method hinges upon the amount of overlap between the volumes Vq and V: 
the tighter Vq bounds V, the more the hits and hence the better the quality of the estimate. Conversely, bounding 
volumes that have little overlap with V give rise to more misses than hits and hence to large errors in the estimate. 
This poses a major obstacle to the implementation of such methods, especially in cases where it is difficult to guess 
the precise location of the volume and its boundaries, often leading to the design of an unnecessarily conservative 
(large) reference region, and hence to very inefficient estimates of V. 

Despite its simplicity, the above example captures the central issues pertinent to integral estimation problems in 
general, for which more sophisticated methods exist [TH3]- Notably, in most efficient Monte Carlo techniques for 
normalizing constant or free energy estimation, a suitable family of intermediate integrands is required to interpolate 
between the two desired integrands (or between the available reference system and the integrand of interest, in absolute 
integral estimation) , the extent of overlap between the integrands dictating the efficiency of the methods much like in 
the example above [3HS] ■ Such requirements become particularly daunting when the integrands in question are highly 
concentrated in unknown and disparate regions of configuration space, as is typically the case in most interesting 
physical problems, thereby demanding substantial user insight prior to the applicability of such methods. 

In the present contribution, a new family of MC integration strategies that can greatly alleviate the demand for 
such types of insights will be introduced. As expressed by the central results underpinning these methods, Eqs. Q 
and ([7]), the integrals of interest [Z, Eq. 0) are computed individually, thus fundamentally departing from the 
aforementioned alternatives that rely on ratios of similar integrals. At the core of these new strategies is the use of 
integration spaces of enlarged dimensionality ( "replicas" ) , a concept already widely invoked to speed up Markov chain 
simulations (parallel tempering [7j, evolutionary Monte Carlo pQ, etc.), combined with the replacement of reference 
integrals by normalized transition matrices [2] (also known as transition functions [T] or kernels [5] in the Markov 
chain literature). Two main variants will be presented: the first, more physically intuitive, requires a fluctuating 
number of replicas (Fig. T] bottom panel); the second, more abstract but also more easily adaptable to existing replica 
simulations, uses a fixed number of replicas along with "virtual" insertions/deletions of replicas. These two variants 
are complementary to each other much like grand-canonical Monte Carlo (GCMC) and Widom's test particle insertion 
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method in simulations involving chemical potentials [2J. 



II. THEORY 



To see how in principle the simultaneous use of multiple configurations (replicas) allows for the absolute computation 
of integrals, let us go back to the volume estimation problem of Fig. [T] As illustrated in the bottom panel of that figure, 
the volume V of interest can be estimated by equilibrating it with an infinite reservoir of ideal gas particles (replicas) 
of density po, and measuring the average number (N) of particles inside V, to obtain V = (N)/pq. This physical 
idea can be implemented computationally by means of a generalized form of the usual grand-canonical Monte Carlo 
method [2J, where the particle reservoir is at density po (or, equivalently, at the corresponding chemical potential), 
and attempted particle insertions/deletions take place in the neighborhood of an existing particle rather than inside 
a fixed volume that bounds V, that neighborhood being defined by a transition matrix T{x'\x). These are the central 
ideas that motivated the development of the two versions described below. 

In general, we would like to estimate integrals of the type 

Z = j dx 7r(x), (1) 
Jn 

where is the support of the integral, and tt(x) is a positive-definite function (see below, however) of the (i-dimensional 
vector x; e.g. ir(x) — e^ E ^ for most physical problems with energy function E(x). Partition functions of discrete 
systems, i.e. Z = ^2 x€ fi^(x), can be dealt with in an analogous fashion. We are given a (normalized) transition 
matrix T(x'\x), such as those routinely used in the trial part of the Metropolis algorithm [2J; for example, T(x'\x) 
could be a uniform probability distribution up to some distance \x' — x\ from x (see dashed circle in Fig. Ill), or a 
Gaussian function centered about x. As with Metropolis and other Markov chain simulations, the width of these 
distributions can be chosen during a preliminary run of the simulation (see below) . 
For non-positive definite integrands, one can invoke the identity 

Z = Z M <sgn(7r(a;)))| 7r |, (2) 

where Z is defined by Eq. Z\^\ = J n dx \tt(x)\, and the average of the sign function of tt(x) is with respect to points 
x sampled from |tt(:e)|. Provided they exist, both quantities on the right hand side of this identity are immediately 
available from the methods below. 



A. Varying number of replicas 



For the first version of the method, we would like to simulate a system of replicas in contact with a reservoir of ideal 
(non-interacting) replicas of specified density p , such that each replica x^ in 51 independently samples the distribution 
n(xi). The corresponding grand-canonical partition function is thus 
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Note that, unlike the traditional case where < N < oo, here the sum starts at N = 1, as in the present method (see 
below) replica insertions/deletions take place in the neighborhood of at least one existing replica. Provided we can 
simulate according to this partition function, the desired integral Z can be found from the equality 

y=i+^, (4) 

which follows straightforwardly from Eq. ^ by computing each moment of separately. (Alternatively, one can 
use (N) = pqZ/(1 — e~ PoZ ) or any other relationship between moments of N and Z, and numerically solve the 
transcendental equation for Z; the question of which Z estimator is more efficient will be left for future studies). 

An algorithm that samples according to the above grand-canonical partition function goes as follows. At the 
beginning of the algorithm we are given at least one point x\ that belongs to f2; let us assume the general case where 
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we already have JV replicas in f2, and let us denote the vector of coordinates of the replicas by x N = (x±, . . . , Xn). We 
then decide whether to insert or remove a replica, typically — but not necessarily (see Methods section) — with equal 
probability. If the decision was an insertion, we sample a new replica coordinate xat+i from T(xpj + i\xi), where Xi is 
a randomly chosen coordinate from the existing JV, and accept its insertion with probability 

v (x N+1 \x^-mm\l N Po< x n+i) ] (S) 

Pacc( 1 ^ ~ | ' jv + 1 y,*Li T{x N+1 \xi) J ' (5) 

except when xm+i lies outside f2, in which case an immediate rejection takes place. Similarly, if the decision was to 
try a deletion, we randomly pick a replica Xj from the existing JV, and accept its deletion with probability 

P^ix™- 1 ^) = min \ 1, ^ ; \ " } , (6) 
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where x N ~~ x is x N excluding Xj, and the sum over i excludes i = j. An exception is the case where only one replica 
remains, which is always rejected. A proof that this algorithm samples according to Eq. Q follows by detailed 
balance (see Methods section). When the replicas correspond to particles inserted uniformly in a fixed volume V, i.e. 
T(x'\x) = l/V, the above acceptance probabilities reduce to those of Ref. [1| (with the due mappings between p and 
chemical potential, and between tt(x) and the Boltzmann factor). Of course, it is also possible to perform ordinary 
7r-preserving Monte Carlo moves on each replica before attempted insertion/deletions [5] (Fig. [3] uses this idea). 

By repeating the above procedure a number of times, the simulation will eventually equilibrate, and the number of 
replicas will fluctuate about its mean value (JV). If the equilibration is too slow, i.e. too few replica insertions/deletions 
are accepted, the width of the distribution T{x'\x) about x can be adjusted accordingly during a preliminary run, 
in a fashion analogous to what is done in Metropolis Monte Carlo to keep the rate of accepted trial moves within a 
reasonable range [2]. Likewise, if the number of replicas starts to grow beyond one's computational capabilities, or 
diminish until it hardly departs from unity, po can be adjusted so that (JV) is a reasonable number consistent with 
one's computing power. In practice, for many-particle systems with extensive free energies (i.e. Z ~ e ), it is best to 
write po = e^ ^ and adjust ^ < instead. 

Note that the present algorithm generalizes standard grand-canonical simulation [2J in two crucial ways. First, it 
inserts and removes entire replicas of the system of interest as opposed to individual particles of a many-body system. 
This conceptual difference is essentially what allows one to relate moments of JV to the partition function of the 
system of interest (Eq. Q). Second, the replicas are introduced in the neighborhood of an existing replica instead of 
inside a fixed region, as prescribed by the transition matrix. This allows for efficient simulation when the integrand 
of interest is sharply peaked about unknown regions of configuration space (cf. Fig. [5| . Although the use of arbitrary 
transition matrices in grand-canonical simulations is known in the mathematical literature [S], to our knowledge the 
use of such ideas for integral/partition function estimation is new. 



B. Fixed number of replicas 



The second version of the replica gas method introduces two important advantages. First, the number of replicas 
is constant as opposed to fluctuating, making it more convenient for parallel computing architectures, and second 
the replicas can be simulated at different temperatures. These features also make the method easily implemented in 
existing parallel tempering (replica exchange) simulations, thereby benefiting from the greatly enhanced equilibration 
rates of these simulations [TJ. The integrals of interest can then be estimated by computing two separate averages 
involving both the integrand ir(x) and the transition matrix T(x'\x), as described below. 

To introduce the method in its simplest form, let us first assume that only two replicas exist (JV = 2), each of them 
independently sampling the distributions tt and fr, via e.g. Metropolis. The integral of interest is Z, as in Eq. (0, 
and Z — J n dx7r(x) is an auxiliary integral; the auxiliary distribution tt is arbitrary (for example, it could be n itself), 
but in typical applications it corresponds to n at a higher temperature, i.e. tt (x) = tt (x), where < (3 < 1. Then 
the following identity holds: 

_ J n dx[Tr(x)/Z] J Q dx' T(x'\x) ■ tt(x') 

J n dx[TT{ X )/Z] J Q dx'[TT( X ')/Z] ■ T(X'\X) 

~ (T(x>\x))-> {) 
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where the average (0(x, x'))f >g of an observable 0(x,x') means that configurations x (x') are sampled from the 
distribution / (<?). Thus, the numerator in the above result requires x to be sampled from tt{x) while x' is sampled 
from T(x'\x) ("virtual replica insertion," in analogy with Widom's method 2 ), and for each such pair of configurations 
one computes the value of 7r(a;'). Likewise, for the average in the denominator, x is sampled from tt(x) while x' is 
sampled from tt{x'), and for each pair one evaluates T(x'\x) ("virtual replica deletion"). In the limit of infinite 
samples, the ratio of these two averages converges to Z as expressed by Eq. 

For simulations with multiple replicas at different temperatures (/3 1; . . . , /?jv), one can simply combine (i.e. sum) 
the above result for each pair of replicas. Thus, for the Ising model results in Fig. |4j the equation 

m) = v m J p " (8) 

was used. The sums run over each replica j at temperature /3j, except j — i. The energy function E{x) for a spin 
configuration x is the usual Ising model function E{x) = — Y2(k a XkXi with periodic boundary conditions [9]. The 
transition matrix adopted T(x'\x) generates a new spin configuration x' by flipping each spin of x with probability 
Pfiip. Thus, T(x'\x) — (pflip)"^ ~ K "(1 — Pmp) N ~^ x where || x' — x | is the Hamming distance (number of spin 
mismatches) between the configurations x and x' , and N is the total number of spins. 



III. RESULTS AND DISCUSSION 



For illustrative purposes, the results of a volume estimation problem in two spatial dimensions using the version with 
varying number of replicas are reported in Figure [2] This example was chosen due to its infinite support, a property 
that would render the use of importance sampling methods difficult, as they would require the design of a non-trivial 
reference volume Vo with similar support and known quadrature (recall that in principle we do not know where the 
integrand is concentrated, or where its boundaries are). The present method performs well in such circumstances 
without any prior information concerning the support of the integrand, by using a simple uniform transition matrix 
T{x'\x) (Fig. [2j dashed square). 

As an application to integrals more general than simple volumes, in Fig. [3] a representative estimate of Z = 
dxsin(x)/x is shown. Note that this integrand is non-positive definite, so Eq. (|2|) was used. The quantities on 
the right hand side of that equation were estimated using the varying number of replicas version of the method, with 
the positive-definite integrand 1 7r(aj) j = | sin(a;)/a;|. 

In Figure [4] the partition function of the two-dimensional Ising model is computed to demonstrate the version with 
fixed number of replicas (similar results are obtained with the non-fixed version). At low temperatures, i.e. ordered 
states, the replicas are densely localized about the spin-up and spin-down states, and hence a local transition matrix 
T(x'\x) is sufficient to ensure efficient convergence of the averages. Conversely, for higher temperatures close to the 
disordered state and above, the relevant configuration space — and hence the spread of the replicas — grows beyond 
the reach of the local transition matrix adopted, thus causing the averages to converge more slowly (cf. right side of 
Fig.§. 

To understand such convergence issues in more detail, consider for simplicity Eq. (JT]) when tt = n (see below for the 
version with varying number of replicas). In order for the averages in Eq. ^ to converge efficiently, the transition 
matrix T(x'\x) has to be such that: 

(a) Most configurations x' sampled from T(x'\x) fall in typical regions of tt for any typical configuration x sampled 
from 7r (so that the numerator is not dominated by those rare configurations with high values of tt(x')); 

(b) Most configurations x,x' sampled from 7r fall in typical regions of T(x'\x) (so that the denominator is not 
dominated by those rare events that cause T(x'\x) to be of appreciable value). 

In the Ising model example of Fig.[4j where T(x'\x) typically flips only a few spins of x, the low temperature estimates 
converge faster as typical spin configurations only differ by a few spins, thereby satisfying both requirements, while at 
higher temperatures close to T c and above, any two typical spin configurations differ by a substantial number of spins, 
and hence the requirement (b) is violated. Adding more replicas, increasing the value of pm p for higher temperatures, 
or using non-local transition matrices (such as those of cluster algorithms [5]) can alleviate the problem, but such 
issues will be left for future studies. 

Note that analogous convergence issues arise in the version with varying number of replicas. Indeed, as can be seen 
by inspection of Eqs. ^ and ([6]), the acceptance probabilities for insertion and deletion are affected by the choice 
of T(x'\x) much like the averages in Eq. Q are affected by criteria (a) and (b) above. (A separate issue is how well 
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the replicas explore the energy landscape. In the fixed number of replicas version, different temperatures are used 
to overcome energy barriers. Although replica exchange operations can be combined with the varying number of 
replicas method, as a proof of concept for the convergence issues above, it suffices to start a population of replicas 
that populate spin- up and spin-down states equally). 

As illustrated by the above example, the most attractive use of the present methods lies in problems where the 
integrand is sufficiently localized, so that a general-purpose, local transition matrix can be used to yield efficient results 
with moderate numbers of replicas. These are precisely the problems for which existing importance sampling-based 
methods are most demanding, and thus these methods can be seen as complementary to each other (see Fig. [5]) . It 
should be noted that the so-called "flat histogram" Monte Carlo methods [TU1 E] are also able to bypass some of the 
difficulties with importance sampling strategies in some cases, especially for discrete systems. However, the required 
human input and scope of such integration methods are rather different: they require the existence and knowledge 
of suitable order parameters and their ranges (this being particularly difficult for entropic problems, such as that 
of Fig. [2j, knowledge of ground state degeneracies, and for continuum systems suffer from systematic errors due to 
discretization schemes, although attempts to alleviate some of these problems have been put forward ;.12j. 

In summary, the present contribution has introduced two variants of a novel Monte Carlo strategy for estimating 
integrals and partition functions. Both versions can be seen as complementary to existing importance sampling or free 
energy methods [U [6] , in that their utility is generally best when the integrands are concentrated in relatively small 
and unknown regions of configuration space. Both continuum and discrete systems are equally amenable to their use. 
By shifting focus from importance sampling functions to transition matrices, it is expected that these methods will 
encourage a change of paradigm in Monte Carlo integration. 



IV. METHODS 



In this section it will be shown that Eqs. ^ and Q satisfy the detailed balance condition 

P (X N ) ■ Ptr (x N+1 \x N ) - Pacc (x N+1 \x N ) = P(X N+1 ) -Ptr(x N \x N+1 ) - Pacc (x N \x N+1 )- (9) 

According to the grand-canonical partition function Eq. (p}]), the probability of observing the microstate x N is given 
by 

p{x N )(x P0<^ N ^M ^ (1Q) 

where the proportionality constant does not depend on the replica coordinates or N. Note carefully the difference 
between the distribution of the labeled microstate x N , corresponding to replica with label "1" being at xi, "2" at 
X2, etc, and that of the unordered set of coordinates {x N } — {xi, . . . ,xjy}, corresponding any replica being at X\, 
another arbitrary replica at x%, etc. This probability is given by 



V{{X N }) = ^P(^) « p»ir( Xl ) ■ ■ ■ tt(x n ), (11) 
p 



where ~^2 P is the sum over all possible permutations of x\, . . . ,xn- Since the replicas are indistinguishable, we have 
p({x }) = Nlp(x ), as in above. Of course, it is possible to use either description (labeled or unlabeled), provided 
the correct probability distributions are used (Eq. (10) or Eq. (11), respectively). In this section, following [T3], we 
will only show the proof using labeled states, i.e. Eq. (10). It is a simple exercise to modify the development below 
for the unlabeled case; the acceptance probabilities, of course, are unchanged. 

Our acceptance probabilities are of the Metropolis-Hastings type, which by construction satisfy detailed balance. 
In the present notation, the formulas are 

Pacc(£ + \x )=mm<M, — r n\ ( ' ( 12 ) 

[ Ptr{X N+l -\x N ) p(x N ) J 

and analogously for p acc (x N \x N+1 ). According to the insertion/deletion algorithm described before Eq. ^6|, the trial 
probabilities for going between the states x N = (x\, . . . , xn) and x N+1 = (xt, . . . , xn, £) are given by 

(xi,...,x N ) s s (x 1: . . . ,x N ,£), (13) 

(i-?)-7vtt 



() 



where pt r (x N+1 \x N ) is given by the expression above the arrows, and pt r (x N \x N+1 ) by the one below them. In the 
insertion trial probability, q = 1/2 is the probability to try an insertion as opposed to a deletion, 1/(N + 1) is the 
probability that the new coordinate £ is inserted in a given slot of the vector x N+1 (in the above case, the last 
slot), 1/N is the probability to pick coordinate Xi as reference, and T(£\xi) is the probability to sample the candidate 
position £ given the chosen reference. In the deletion trial probability, (1 — (?) = 1/2 is the probability to try a deletion, 
and l/(iV + 1) is the probability that the replica at £ will be chosen for attempted removal among the existing ones. 
Plugging these trial probabilities in Eq. (12), we obtain Eq. QSJ) (an analogous procedure gives Eq. (J6j ) - The case 
where q ^= 1/2 can be easily taken care of by modifying the acceptance probabilities accordingly. 

Alternatively, detailed balance can be directly proven by plugging the trial probabilities in Eq. ( 13 ) and the accep- 
tance probabilities given by Eqs. (5) and ^ into Eq. 
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FIG. 1: Monte Carlo estimation of a volume V (shaded region) by means of sampling from a reference volume Vo (top), and 
equilibration with a hypothetical, infinite reservoir of ideal gas particles at density po (bottom). In the former, one draws 
random points uniformly from Vo and counts the fraction / that lands in V, to obtain V = Vof. In the latter, one performs a 
grand-canonical Monte Carlo simulation [2] at reservoir density po, and monitors the average number (iV) of ideal gas particles 
in V; upon equilibration, the density of particles in V equals that of the reservoir, and thus V = (N) / po- (In practice, due to the 
constraint N > 1, this formula for V needs to be modified slightly; see Eq. Q). Each particle corresponds to a point ("replica") 
residing in V, and attempted replica insertions/removals take place in the neighborhood (dashed circle) of an existing replica 
x, defined by the transition matrix T(x'\x) of the method. A version of the algorithm with fixed number of replicas — possibly 
at different temperatures — is also described in the text. 




FIG. 2: Estimation of a two-dimensional volume Z (yellow region) using the replica gas method with varying number of 
replicas, Eq. Q. The volume Z is defined by the region \x% | < e ~ lxil + 1 for |:ti | > 1, and \x 2 \ < 1 - ln|a;i| for < 1. The 
points correspond to the replica configurations at the end of one simulation, and the dashed square defines the boundaries of 
the adopted transition matrix T(x'\x) (uniform distribution, each side of length unity) for the particular configuration x shown. 
Inset: Histogram of 20 independent estimates of Z using the replica gas method with 10 6 attempted insertions/deletions, and 
po — 5. The exact value of Z, obtained by analytic quadrature, is Z = 12. 
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FIG. 3: An illustrative non-positive definite integrand, n(x) — sin(:r)/:r. The scale in the center of the graph corresponds to 
the parameter A = 5 in the uniform transition matrix T(x'\x) = 1/A for \x' — x\< A/2 (zero otherwise) adopted for the results 
shown in the inset. Inset: Running estimate of Z = dxTv(x) using Eq. Q2p, where Z\ v \ is estimated via the replica gas 
method with varying number of replicas, Eq. Q. The mean sign function of ir required by this last equality, (sgn(7r(a;)))| vr |, is 
also obtained from this run, by averaging sgn(7r(a:)) over all replicas x during the simulation. The dashed red line is the exact 
result Z = 7T = 3.141592.... For this example, po = 1, and 10 ordinary Monte Carlo moves per replica are performed between 
every attempted insertion/deletion (GCMC step). The ensuing number of replicas fluctuated about N = 10. 
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FIG. 4: Natural logarithm of the partition function Z of the two-dimensional Ising model with N = 32 x 32 spins, according 
to the replica gas method (circles), and the exact Kaufman formula [5] (dashed curve). The replica gas results were obtained 
using the fixed number of replicas version, Eq. (8j, with 20 replicas at the temperatures corresponding to the data points shown 
(similar results are obtained with varying number of replicas, see text), with error bars indicating the standard deviation of 8 
independent runs. Importance sampling results were obtained using the ideal (non-interacting) spin partition function Z^ = 2^ 



and Zi s i n g/.Z id = (e 



-/3-ElsmirO) 



)id, with 10 independent configurations x sampled from the ideal reference system (increasing 



this number to 10 does not lead to appreciable changes in the results). For disordered states (i.e. fcsT close to or higher than 
the critical temperature ksT c — 2.269), the partition function is no longer dominated by a small fraction of the configuration 
space, and the replica gas method converges more slowly with the adopted (local) transition matrix (see also Fig. [5j. The 
replica exchange simulation took 10 5 MC steps, with an attempted exchange every 100 steps, where each MC step is a simple 
spin flip. For the transition matrix sampling, paip — 1/N (cf. discussion after Eq. (ISJ) ) . 
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FIG. 5: Comparison of the merits of importance sampling (left) and replica gas (right) methods for sparse (top) and localized 
(bottom) integrands. Integrands are represented by their densest regions (shaded curvy shapes). For physical systems, sparse 
integrands correspond to Boltzmann factors at high temperatures, while at low temperatures the integrands tend to be localized 
in a small fraction of configuration space (e.g. magnetized spin systems, crystals, proteins in their native state, etc). In 
importance sampling, one typically has at their disposal a general-purpose sparse reference system (polygon on the left) 
such as an ideal gas, which is generally sufficient to ensure proper sampling at high temperatures, but not at lower ones. 
Conversely, in replica gas methods one typically has at their disposal a local transition matrix (small boxes on the right) that 
is generally sufficient to sparingly "cover" the integrand of interest (cf. efficiency criteria (a) and (b) discussed in the text) at 
low temperatures, but not at higher temperatures. 



